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Abstract. In this work, we investigate a particular class of shape optimization problems under uncertainties 
on the input parameters. More precisely, we are interested in the minimization of the expectation of a 
quadratic objective in a situation where the state function depends linearly on a random input parameter. 

This framework covers important objectives such as tracking-type functionals for elliptic second order partial 
differential equations and the compliance in linear elasticity. We show that the robust objective and its 
gradient are completely and explicitly determined by low-order moments of the random input. We then 
derive a cheap, deterministic algorithm to minimize this objective and present model cases in structural 
optimization. 
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1. Introduction 

Over the last decades, shape optimization has been developed as an efficient method for designing devices 
which are optimized with respect to a given purpose. Many practical problems in engineering lead to 
boundary value problems for an unknown function which needs to be computed to obtain a real quantity of 
interest. For example, in structural mechanics, the equations of linear elasticity are usually considered and 
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solved to compute e.g. the leading mode of a structure or its compliance. Shape optimization is then applied 
to optimize the workpiece under consideration with respect to this objective functional. We refer the reader 
to [2 na Eoi Ea [29] and the references therein for an overview on the topic of shape optimization, which 
falls into the general setting of optimal control of partial differential equations. 

Usually, the input parameters of the model, like the applied loads, the material’s properties (typically the 
value of the Young modulus or Poisson ratio) or the geometry of the involved shapes itself are assumed to 
be perfectly known. However, convenient for the analysis of shape optimization problems, this assumption 
is unrealistic with respect to applications. In practice, a manufactured device achieves its nominal geometry 
only up to a tolerance, the material parameters never match the requirements perfectly and applied forces 
can only be estimated. Therefore, shape optimization under uncertainty is of great practical importance but 
started only recently to be investigated, see e.g. mmiHiisiiiiiiniiisi for related results. 

Two approaches are available in the context of optimization under uncertainty, depending on the actual 
knowledge of the uncertain parameters. On the one hand, if no a priori information is available short of an 
upper bound on their magnitude, one usually considers a worst-case approach. On the other hand, if some 
statistical information on the distribution of the unknown parameters is given, one can study the statistics of 
the objective, which depends on the random parameters through the state equation. Notice that in this case 
the state function is a random field and so the objective itself becomes random - the value of the objective 
depends on the design variables and on the random variable. 

One is usually primarily interested in stochastic quantities of the objective such as its expectation. When 
this crude average is not sufficient, one may consider a weighted combination of the expectation and the 
standard deviation in order to limit the dispersion of the objective values around its expectation. Finally, 
one sometimes also considers the probability that the objective exceeds a given threshold. This last objective 
usually stands for constraints. 

In the present article, we address the following problem: given a partial statistical description of the 
random loading, design an efficient algorithm to minimize the expectation of the objective. 

We restrict ourselves to a special class of problems, namely those involving a quadratic shape functional 
for the state function which is defined by a state equation with a random right-hand side. This in particular 
means that the random state depends linearly on the random input parameter. Our theory covers important 
shape functionals like the Dirichlet energy and quadratic tracking-type functionals in the context of the 
Laplace operator. Besides, the compliance functional in linear elasticity belongs to the important members 
of the class of functionals under consideration. 

Our main message is the following: for objectives of the class under consideration, whose expectation is to 
be minimized, all quantities for performing a gradient-based shape optimization algorithm can be expressed 
deterministically, i.e., without any approximation. We only need access to the random parameter’s first and 
second moment. The main object is the two-point correlation function Cor(it) of the state function u. It is 
the solution of a tensor-product type boundary value problem with the random right-hand side’s two-point 
correlation as right hand side. As a consequence, both, the expectation of shape functional, as well as the 
related shape gradient can explicitly be determined and efficiently be computed just from the knowledge 
of the random right hand side’s expectation and two-point correlation function. This fact is of tremendous 
importance for applications: it is completely unrealistic to have access to the law of the random loadings, 
whereas the knowledge of its expectation and of its two-point correlation function seem to be much more 
reasonable assumptions. We therefore end up with a fully deterministic, efficient algorithm of similar cost 
as for classical shape optimization when no uncertainties are taken into account. 

This article is organized as follows. First, we present in Sectionj^the leading idea to reduce the stochastic 
problem to a deterministic one, relying on a very simplified model in finite dimension. We introduce the 
tensor formulation that is the keystone of the subsequent calculations. Then, in Section we present the 
shape calculus which we shall use and adapt the idea of Sectionj^to this more complex setting. In particular, 
we recall definitions and properties of tensor products on Hilbert spaces. We then apply in Section the 
obtained method to three significant examples in the context of the Laplace operator and the equations of 
linear elasticity. Finally, we explain in Section how to design efficient numerical methods to solve the 
corresponding optimization problem. The main point is to remark that the numerical resolution of the high¬ 
dimensional boundary value problem which defines Cor(rt) can be avoided if desired. An appropriate low-rank 
approximation scheme allows to reduce this computation to the resolution of some classical boundary values 
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that can be solved thanks to a standard toolbox. This leads to a non-intrusive implementation of the 
proposed method. We conclude in Section]^ with numerical examples concerning the robust optimization of 
the compliance of a mechanical structure. 


2. Formal presentation of the main idea 

In this section, we formally outline the main idea of our approach in a finite-dimensional setting where 
calculations can be performed in an intuitive way by using only elementary algebra. To this end, let "H be a 
vector space of designs h, whose performances are evaluated by a cost function C{h,uj) which depends on h 
via the solution u{h,uj) = of the fV-dimensional system 

( 2 . 1 ) A{h)u{h,ui) = 

In this formula, A{h) S is an invertible matrix of dimension N x N, f{u;) is a (random) vector in K^, 
and w € is an event, belonging to a complete probability space (fl, S,P). The cost function C is assumed 
to be quadratic, i.e., of the form 

C(h, uj) = {Bu{h, oj), u{h, uj)) = B : (u(h, ui) 0 u{h, to )), 

where B € is independent of the design for the sake of simplicity. In this formula, the tensor product 
u 0 w of two vectors v,w £ is the (TV x TV)-matrix with entries {v 0 w)ij = ViWj, i,j = l,...,N, and : 
stands for the Frobenius inner product over matrices. 

The objective function of interest is the mean value of the cost C{h,uj): 

(2.2) M{h)= [ C{h,uj)F{duj) = B : Coiiu){h). 

JQ 

Here, Cor(u)(/i) is the N x N correlation matna0of u{h,uj), whose entries read 

CoT{u){h)ij = / Ui{h,oj)uj{h,uj)F{doj), i,j = l,...,N. 

JQ 

This matrix can be calculated as the solution to the A^^-dimensional system 

(2.3) (-^(^) ® -^(^)) Cor(u)(/i) = Cor(/). 

At this point, let us recall that A{h) 0 A{h) : —)■ is the unique linear mapping such that 

Vu, V S R^, {A{h) 0 A{h)){u 0 u) = {A{h)u) 0 {A{h)v). 

Let us no w calcul ate the gradient of Ai{h). Denoting with ' the differentiation with respect to h, we 
differentiate (2.1 2.2) in an arbitrary direction h to obtain 

(2.4) M'{h)(h) = 2 [ {Bu\h,uj)(h),u{h,oj))F{doj), 

Jn 

where 

A{h)u {h,uj){h) = —A'{h){h)u{h,uj). 

Introducing the adjoint state p{h,uj), which is the solution to the (Wdimensional) system 

(2.5) A{h)^p{h, Lo) = —2B^u{h,u!), 
we derive successively 

2{Bu\h,uj){h),u{h,uj)) = —{u'{h,uj){h),A{h)'^p{h,oj)) = {A' {h){h)u{h,uj),p{h,oj)). 

Hence, we arrive at 

M'{h){h) = ^A'(/i)(h) 0 Cor{u,p){h). 

In this last formula, the correlation matrix Cor(u,p)(/i) can be calculated as the solution to an iV^- 
dimensional system; indeed, using ( 2.1|2.5 ), one has for any event uj that 

{A{h) 0 A{h)'^) {u{h, oj) 0 p{h, uj)) = — (A(h) 0 B) {u{h, uj) 0 u[h, uj)), 


^Unfortunately, the literature is not consistent as far as the notions of covariance and correlation are concerned. In this 
article, the correlation between two square-integrable random variables X and Y is defined as E(Xy). 
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whence the following system for Cor('u,p): 

(2.6) (g) A{h)'^') Cor{u,p){h) = — {A{h) (g B) Coi{u){h). 

These considerations show that both, the objective function and its gradient, can be exactly calculated 

from the sole datum of the correlation matrix of / (and not of its law!). 

At this point, one may wonder about the practical interest of the foregoing computations since the 
systems (2. ^ seem difficult to solve (see however [27)1. The main idea consists in performing a low-rank 

approximation of the correlation matrix Cor(/): 

m 

Cor(/) « ^/i(g/,, m<^N. 

2=1 


Then, formula (2.3) leads to the calculation of a convenient approximation of Cor(n)(h) in accordance with 

m 

Cor(u)(h) « ^ Ui{h) (g Ui{h), 


i=l 


where Ui{h) arises as the solution of the system 

(2.7) A{h)u^{h) = fi. 


Similarly, by (2.61, one has 


Cor(M, p){h) » Ui{h) (g p^{h) 


i=l 


with 


( 2 . 8 ) 


A{hfp^{h) = -B^u,{h). 


Hence, calculating A4{h) and its derivative A4'{h) amounts to solving (only) m systems of the form (2.7) 


and m systems of the form (2.8). 


Remark 1. 


Cost functionals C of the design involving a linear term of the form i{u(h,uj)) can be considered in 
the same way (see Section [4.1.2 ). The corresponding mean value also involves the mean value of u, 
E(u)(/i) := \^ u{h ,oj) P((icLi). 

Formulae ( |2.2|2.4 ) show explicit expressions of A4{h) and M'{h) only in terms of the correlation 
Cor(/), which is quite appealing for at least two reasons. First, Cor(/) may be imperfectly known (in 
realistic applications, it is often reconstructed from observations by statistical methods). Second, as 
we have just discussed, it is often desirable to approximate it, so to ease numerical computations. In 
either situation, these formulae allow to measure directly the impact of an approximation of Cor(/) 
on M{h) and 

An alternative approach for the calculation of A4(h) and M'(h) consists in computing a truncated 
Karhunen-Loeve expansion of /(cu), i.e., /(w) « with {fi} being orth ogon al vectors 

in and {^i} being uncorrelated random variables. Injecting this expression in (2.1) yields an 
approximation of u in accordance with u{h,uj) « where the Ui(h) are given by 

(2.7). Then, using the quadratic structure of the cost C allows to conveniently approximate M.{h) 
and A4'{h), leading to similar formulae. Doing so is however less efficient than the proposed approach 
for at least two reasons. On the one hand, calculating the Karhunen-Loeve expansion of a random 
field is rather involved in terms of computational cost. In contrast, the proposed method in this article 
relies on any low-rank approximation of the correlation Cor(/). On the other hand, estimating the 
error entailed by such a process is awkward, since it does not rely on the direct connection between 
M{h), M'{h) and Cor(u), Coi{u,p) (it passes through the approximation of u{h,uj) itself, which is 
not directly involved in their expressions). 
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3. Shape optimization setting 


Extending the framework presented in the previous section to the infinite dimensional setting, and more 
specifically to that of shape optimization, demands adequate extensions of the notions of random variables, 
correlation matrices, etc. At the centre of these generalizations lies the notion of tensor product between 
Hilbert spaces, about which this section gathers some useful material for the reader’s convenience. 

3.1. Differentiation with respect to the domain: Hadamard’s method. 

Several notions of differentation with respect to the domain are available in the literature. In this article, we 
rely on Hadamard’s boundary variation method (see e.g. [3 nil US]), which advocates to describe variations 
of a bounded, Lipschitz domain C of the form (see Figure]^: 

(3.1) De = {I + e){D), 1. 



Do 


Figure 1. One variation Dg of a shape D in the framework of Hadamard’s method. 


Definition 1. A function J{D) of the domain is said to be shape differentiable at D if the underlying 
functional 9 i-A J{Dg) which maps K"’*) into K is Frechet-differentiable at 9 = 0. The shape 

derivative 9 i-A J'{D){9) of J at D is the corresponding Frechet derivative, so that the following asymptotic 
expansion holds in the vicinity of 9 = 0: 

J{Dg) = J{D) + J'{D){9) + o{9), where ,,,,, ^ 0. 


In practice, shape optimization problems are defined only over a set lAad of admissible shapes (which 
e.g. satisfy volume or regularity constraints). To ensure that variations (3.1) of admissible shapes remain 
admissible, one usually imposes that the deformations 9 he in a subset Qad C K"^) of admissible 

deformations. 

In the following, we implicitly and systematically assume that the sets lAad and Qad contain shapes or 
deformations with sufficient regularity to legitimate the use of the classical formulae for the shape derivatives 
of the considered functionals. We refer e.g. to m for precise statements on these issues. 


3.2. Tensor products and Hilbert spaces. 


In this subsection, we collect some dehnitions and properties around the notion of tensor product of Hilbert 
spaces. A more thorough exposition can be found in [BUS]. In what follows, we consistently assume all 
the Hilbert spaces to be separable. This is merely a matter of convenience and most of the forthcoming 
dehnitions and results hold also in the general context. 
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Definition 2. Let {Hi, {•)hi), {H 2 , {■)h 2 ) be two (separable) Hilbert spaces. Then, for any hi G Hi, /12 G H 2 , 
the pure tensor hi 0 /12 is the bilinear form acting on Hi x H 2 as 

(3.2) V((/5i,</ 52) G i?i X i?2, {hi®h2){(pi,y}2) = {hi,^Pi)Hi{h2,T2)H2- 
The vector space of all pure tensors 

= span{/ii (g)/i 2 , /ii G i?i, ft .2 G i? 2 } 

has a well-defined inner product H x H —>■ K which is determined by its action on pure tensors 

(3.3) Vtpi ® (^2, '01 0 ^^2 G (</9l g) (/52, 01 0 02) = (‘/?l,0l)jli(v^2, 02 ) 1/2 

and extended to TL by bilinearity. In particular, this inner product does not depend on the choiee of the 
decomposition of elements of TL as finite sums of pure tensors used to calculate it. The tensor product 
Hi g) H 2 is finally the Hilbert space which is defined as the completion of TL for (•, •). 


Remark 2. Alternative definitions can be found in the literature, e.g. relying on the notion of Hilbert- 
Schmidt operators. 


Proposition 1. Let Hi, H 2 be Hilbert spaces and let {0i}jg|!^, {0j}jgis} be associated orthonormal bases. 

(1) The set {(pi g) 0^}^ is an orthonormal basis of Hi g H 2 for the inner product defined in 

(2) For any h & Hi® H 2 , there exists a unique family of elements in Hi such that 

00 

h = ^ Ui ® 0i. 
i=l 

The main purpose of tensor algebra is to transform multilinear expressions into linear ones. Before stating 
a version of the universal property of tensor products in the Hilbert space context, let us recall that a weak 
Hilbert-Schmidt mapping b : Hi x H 2 — >■ K between Hilbert spaces Hi, H 2 , {K, Qk) is a bilinear mapping 
with the property that there is a constant c > 0 such that 

00 

VuGAT, ^ 1(6(0*, 0j),M)icP < c||u||^ 

for given (then any) orthonormal bases { 0 i}jgf^ of Hi and { 0 ^}^.^^^ of H 2 , respectively. 


’3.3). 


Proposition 2. Let Hi, H 2 , K he Hilbert spaces, and b : Hi x H 2 K be a weak Hilbert-Schmidt mapping. 
Then, there exists a unique bounded operator! : Hi®H 2 —t K such that the following diagram is commutative 


Hi X H2 



where the mapping Hi x H 2 — t Hi g H 2 is simply {hi, 62 ) 1 — >■ 61 g 62 . 


We now come to the following very important identification of L^ spaces taking values in a Hilbert space. 

Proposition 3. Let {LI, p) be a measured space, and H be a Hilbert space. Then, the mapping 

L^{Vt, p)xH 3 {f, h)^fhG L^{Vt, p, H) 
induces a natural isomorphism L^{n,p) g R ~ Lf{Ll,p,H) between Hilbert spaces. 
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Example 1. In the particular case that H = Lp'{D, v), where {D, v) is another measured space, Proposition 
[^yields the isomorphism /i) 0 L^{D, v) « x D, ly), where ly stands for the usual product 

measure of /r and v on x D and the above identification is supplied by 

Vu S fi), V G L^{D, v), (u ® v){x, y) = u{x)v{y), x G ft, y G D. 

In the following, we consistently employ this identification. 

3.3. First and second moments analysis. 

In this section, we slip into the probabilistic context, so to speak, relying on the framework of [miMi. Hence, 
let (H, E,P) be a complete probability space and let H he a (separable) Hilbert space. To keep notation 
simple, we omit to mention the measure P on H when the context is clear. 

Definition 3. (1) The mean value E : ® H ^ H is the unique linear and bounded operator which 

satisfies 

G uGH, E(e <^u)= (^J u. 

(2) The correlation Cor(u,u) G H ® H between two elements u,v G ® H is defined as 

OO 

Cor(u, v) = Ui ® Vi, 

i=l 

where u = ® ® decompositions of u and v supplied by Propo¬ 
sition^ according to an orthonormal basis of Cor(u,v) is moreover independent of 

the basis used to perform the above construction. 

(3) The function Cor{u,u) is simply denoted as Cor(u) and called the (two-point) correlation of u. 

This terminology is consistent with the usual definitions of the mean and correlation of random fields. 
Indeed, if ZI C is a domain (equipped with the usual Lebesgue measure) and H — Lfi{D), then it is easily 
seen that the mean value E(tt) G L^{D) of an element u G L^(fl) 0 L'^{D) — (ft x D) is 

E(m) = / u(-, w) P(fia;), 

Jn 

and that the correlation Cor(M,u) G L'^{D) 0 L'^{D) cx L'^(D x D) oiu,v G 0 L^{D) is given by 

Coi{u,v){x,y) = / u{x,Ljj)v{y,u})¥{dijj) a.e. (x,y) G D x D. 

Jn 

For further use, we have to give a sense to expressions of the form Cor(u, u)(a;, a;), where u,v G 0 

L‘^{D) and x G D. Note that this is not completely straightforward since Cor{u,v) is only an element in 
Lfi{D X D); hence, it is a priori not defined on null measure subsets oi D x D. Doing so is the purpose of 
the next lemma: 


Lemma 4. (1) Let F 0 Lfi{D x D) be the subspace defined as 

= span |u (g) u, u,v G 

equipped with the nuclear norm 


{ N N \ 

\ \u^\\L^(^D)\\y^\\L^{D): h = '^Ui^Vi, Ui,Vi € {D) \ . 

There is a unique linear and continuous operator 7 : —> L^{D) such that, for any functions 

u,v G T^D), it holds that 

7 (it 0 v)(x) = u(x)v(x) a.e. x G D. 
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(2) Let Ta C L^{D X D) be the subspace defined as 

Tc = span {Cor(u, v), u,v G ® Lfi{D)^ , 

equipped with the norm 

{ N N 

^\\ui\\L^{axD)\\vt\\L2inxD)j h = ^C0I{ui,v^), Ui,Vi € L^{9.)®L^{D) 

There is a unique linear and continuous operator '■ ^ such that, for any functions 

u,v G 'D(fil X D), it holds that 

'yc{CoT{u,v)){x) = Cor{u,v){x, x) a.e. x G D. 

(3) The following diagram is commutative: 


L'^{n X D)'^ L^{D X D) 


la 


Proof. 1. First, note that H-H* does define a norm on J- since, for arbitrary h G J-, one has that ||h|| l 2 (£)xd)< 
||/i||*. What is more, the subspace G := spanjrt® v, u,v G ViD)} is obviously dense in iF for the norm 
11-11*. Define now the mapping j : G by 

N N 

yh = ''^Ui ^ Vi, Ui,ViGT>{D), j{h){x) =''^Ui{x)vi{x) a.e. X G D. 

i=l i=l 

This mapping is obviously well-defined and continuous provided that G is endowed with the norm ll'H*. 
Thus, it is uniquely extended into a linear and bounded operator ■j : T ^ L^[D), which fulfills the desired 
properties. 

2. The proof is completely analogous to that of 1. 

3. The commutation relations obviously hold for smooth functions G 'D{Vl x D), and the general result 

follows by continuity of the mappings at play. □ 


Remark 3. 

• To keep notations simple, we will often omit to mention explicitly the operators 7 and 7 c, and, e.g. 
we write expressions such as (u ® v){x, x) instead of ^(u ® v)(x, x). 

• Analogous definitions and commutation relations, involving different operators (derivatives, etc.), 
hold and can be proved in an identical way when the space H used in the Definition of the 
correlation function is, for instance, H = [L^{D)]'^ or H = H^{D). We do not specify these natural 
relations so to keep the exposition simple. 

4. Applications in concrete situations 

4.1. Quadratic shape functionals in the context of the Poisson equation. 


Let D be a smooth domain and / G L^(R‘^)) be a random source term. For almost any event uj G ft, 

let ud{-,uj) G Ho{D) be the unique solution to the Poisson equation 


(4.1) 


-Au{-,uj) = f{uj) in D, 
u{-,uj) = 0 on dD. 


By the standard Lax-Milgram theory, this equation has a unique solution in Hq{D) for almost all events 
uj G LI, and it is easily seen that ud G L^(L2, Hq(D)). Moreover, owing to the standard elliptic regularity 
theory (see e.g. 0 ), it holds u(-,uj) G H^{D) for almost all events u gLI. 






4.1.1. The Dirichlet energy as cost function. 


The first functional under consideration is the mean value A4{D) of the Dirichlet energy: 

M{D)= f CiD,uj)¥{duj), C{D,uj) = -l f \\VuD{x,uj)\f dx. 

Jn ^ Jd 

The result of interest is as follows. 

Theorem 5. The objective function Ai{D) can he rewritten as 

(4.2) 7\d(£)) = — i f (V 0 V) Cor(zi£))(x, x) da;. 

2 J D 

This functional is shape differentiable over lAad o.nd its shape gradient is given by 

(4.3) '^d e Qad, M'{D){e) =-^ J Cor {ud){x,x) {9 ■ n){x) ds{x). 

In the last two formulae, (V 0 V) : Hq{D) 0 Hq{D) — >• L'^{D) ® L'^{D) and -§^) '■ H'^{D) ® H‘^{D) — >• 

L^{dD) ® L^{dD) stand for the linear forms induced by the respective bilinear mappings 

/ N 7 / \ du dv 

{u, v) Vu • \/v and [u, v) ;—>■ 7;— 7;—. 

on on 

The correlation function Cot{ud) & Hq{D) x Hq{D) ~ Hq{D x D) can he obtained as the unique solution 
of the boundary value problem 

{ —(A (g) A)(Cor(u£))) = Cor(/) in D x D, 

(4.4) < 

y Cot{ud) =0 on d{D x D). 

Proof. We find 

A4(Zl) = -i f f ||VMD(x,a;)|p dxP(da;) 

2 Jn Jd 

=-\ [ /" 7((V (g) V)(u_d(-,w) (g) MD(-,a;)))(x,x) (ixP(dw). 

2 Jn Jd 

Now using Lemma(see also Remark]^, we obtain 

A4(L>) = -i y 7c (V (g V)(u_d(-,w) 0 MD(-,a;))(x,x) P(da;)^ dx 

= y 7c 0 V) ^y ud 0 ud P(dw)^^ (x,x) dx, 


which implies the desired expression (4.21. 


To prove (4.3), we follow the same analysis, starting from the classical formula for the shape derivative. 


for a given event w G fl (see e.g. my- 

C'{D,uj){0) = - 


f 

9 u d (-, w ) 

JdD 

dn 


9 ■ n ds. 


Finally, that ud is the unique solution to the system (4.4) follows from tensorizing the state equation 
(4.1) and applying the usual Lax-Milgram theory (again, see or details). □ 


4.1.2. L'^-tracking type cost function. 

Still in the setting of the Poisson equation outlined above, we are now interested in the L^-tracking type 
cost function 


C{D,uj) = ^J \uDix,uj) - uo{x)\‘^ dx, 
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where u_d is the solution to (4.1), B d D is a fixed subset of D, and uq € is a prescribed function. 

Again, we aim at minimizing the mean value of the cost C, that is 

M{D)= [ C{D,uj)F{dio). 

JQ 

The main result in the present context is the following theorem: 

Theorem 6. The functional M.[D) defined above can be rewritten as 

(4.5) M{D) = ^ I (^Cor{uD){x,x) — 2uo{x)E{u£)){x) + ul{x)) dx. 

2 Jd 

It is shape differentiable at any shape D G Uad with shape derivative given by 

(4.6) '^e G Bad, M'{D){0) = - J Cot{pd,ud){x,x){ 9 ■ n){x) ds{x). 

In this formula, the adjoint state pu G 0 IIq{D) satisfies the boundary value problem 

-Ap(-,a;) =-xb(md(-, w) - uo) in D, 
p{-,oj) = 0 on dD, 


(4.7) 


a.e. uj G fl, 


where xb stands for the characteristic function of B. Moreover, the mean value E(u£)) G IIq{D) which 
enters El is the unique solution of 


(4.8) 


-AE(w) = E(/) in D, 
E(m) = 0 on dD, 


the two-point correlation Cot:{ud) is the solution of (j.j), and the correlation function Cot{pd,ud) G 
Hq{D) (g) IIq{D) can be calculated by solving the boundary value problem 


(4.9) 


— (A g /)Cor(p, u) = —{xb <8) /)(Cor(u) — ug 0 E(m)) in D x D, 

Cor(p, u) =0 on d{D x D). 


Proof. The proof is essentially identical to that of Theorem]^ and hence is not repeated here. Nevertheless, 
let us just point out that (|4.9|) stems from the computation 


-{A0 I)CoT{pD,UD){x,y) = / -ApD{x,u!)uDiy,uj)F{du}) 

Jn 

=-Xb{x) / {ud{x,uj) - uo{x))uD{y,u})F{duj) 
Jn 


= -Xb{x) (Cor(u)(x,y) - uo(x)E(un)(y)) 


for almost all (x, y) G D x D. 


□ 


4.2. Quadratic functionals in the context of linear elasticity. 

We now slip into the context of the linear elasticity system. The shapes D under consideration are 

filled with a linear elastic material with Hooke’s law A given by 

Ve e 5(11^^), Ae = 2//e + A tre/, 

where the Lame coefficients X and ji satisfy y > 0 and A + 2/i/(i > 0. 

The admissible shapes D G Uad are clamped on a fixed subset T^i of their boundaries, surface loads are 
applied on another fixed, disjoint part Tjv C dD, so that only the free boundary T := dD \ (T^i U Tjv) is 
subject to optimization. Accordingly, we shall assume that all the deformation fields 0 G &ad vanish on 
Td UTat. Omitting body forces for simplicity, the displacement ud of D belongs to the space [H^ {D)]‘^, 
where 

H^^{D) = {uG H^{D), s.t. u = 0 on Tij} , 
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and is the unique solution in this space to the boundary value problem 


(4.10) 


—div(Ae(u)) =0 in 

u = 0 on r^), 

< 

Ae{u)n = g on Fat, 
Ae(u)n = 0 on F, 


where e{u) = (Vu + Vu^)/2 stands for the linearized strain tensor. 
The cost function at stake is the compliance of shapes 


C(Acc) 


Ae(uD){x,uj) : e{uD){x,uj) dx = 


g{x,ui) ■ ud(x,uS) ds{x), 


jd Jd 

and we still aim at optimizing its mean value A4(D) = J^C{D,uj) P(da;). Arguing as in the previous sub¬ 
section, we obtain the following result. 


Theorem 7. The above functional Ai{D) can be rewritten in accordance with 


M{D) = 


ID 


{{Acx '■ ey)Cor{u))(x,x) dx, 


where {Ae^ '■ Cy) : [Hf^[D)Y ® [Hy^{D)Y L'^{D) ® L^{D) is the linear operator induced by Proposition 
^from the bilinear mapping 

{u, v) !-)■ Ae{u) : e{v). 

This functional is differentiable at any shape D S Uad und its derivative reads 


V0G0ad, M'{D){e) = - j^{{A 

Here, the two-point correlation function Cor(M) G 
following boundary value problem: 

' (divj; 0 d\\y){Aex ® Acy) Cor(M) = 0 

Cor(u) = 0 

(divj, ® Iy){Aex ® ly) Cor(u) = 0 
{Ix ® divy){Ix ® Acy) Cor(M) = 0 
{Acx ® Acy) Cor{u){nx ® Uy) = Cor((?) 
(diva; (g) Iy){Aex ® Acy) Coy{u){Ix ®ny) =0 
{Ix ® divy)(Aex ® Acy) Coi{u){nx ® ly) =0 
{Acx ® Acy) Cor{u){nx ® Uy) = 0 
{Acx ® ly) Coi{ u){nx ® ly) =0 

{Ix ® ACy) COl{u){Ix ® Uy) = 0 


Cx '■ ey)Coi{u)){x, x){9 • n){x) ds{x). 

\Hf^{D)Y ® \Hyo^D)Y is the unique solution to the 

in D X D, 
on Fd X F/j, 
on D xT £), 
on Fd X D, 
on Fat x Fat, 
on X (Fat U F), 
on (FatUFat) X D, 

on ((Fat U F) x (Fjv U F)) \ (Fjv x Fat), 
on (Fat x F) x Fa), 
on Fa) X (Fat X F). 


Remark 4. All the involved mappings in the foregoing expressions are naturally produced by Proposition 
and we do not make the underlying functional spaces explicit. The subscripts x and y refer to operators 
acting respectively on the first and second component of a pure tensor in a tensor product space. 


5. Numerical realization 

In this section, we now focus on how the previous formulae for the objective functions of interest and 
their derivatives pave the way to efficient calculations in numerical practice. 
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5.1. Computing second moments. 


Without loss of generality, we focus the discussion on the setting of the Poisson equation, as discussed in 
Section 4.1 The expressions ( 4.2|4.3|4.5|4.6 ) involve the mean value E(m) and the correlation Cor(u) of the 
solution u{-,w) to (4.1) and the correlation Cor(u,p) between u and the solution to (|4.7[). 


The quantity E(m) is fairly straightforward to calculate once the mean of the data E(/) is known; indeed, 


it arises as the solution to the boundary value problem (4.8) which can be solved owing to any standard 
finite element method. 


It is however more complicated to compute Cor(u) (or Cot{u,p)) since, in accordance with (4.4), a fairly 
unusual boundary value problem for the tensorized Laplace operator needs to be solved on the product 
domain D x D. This moderately high-dimensional problem can be solved in essentially the same complexity 
as (4.8) if a sparse tensor product discretization is employed as proposed in e.g. [13 [mill]. However, the 


implementation of this approach is highly intrusive, insofar as it demands a dedicated solver. 

One way to get past this difficulty, which is also much simpler to implement, consists in relying on an 
expansion of the two-point correlation function of of the form 


(5.1) 


Cor(/) =^/feO/fc. 


Then, the two-point correlation function Cor(u) can be expressed as 


(5.2) 


Cor( 


= E 


Uk O Ufe, 


where each function Uk is the solution to the Poisson equation ( |4.l| ) with data fk- In other terms, the solu¬ 
tion’s two-point correlation function can be determined from solving (possibly infinite) standard boundary 
value problems. In practice, the expansion (5.1) is truncated so that this process becomes feasible. 


Several possibilities are available when it comes to decomposing Cor(/) as in (5.1). For example, in the 
situation that Cor(/) G L‘^{D x D), the most natural idea is to perform a spectral decomposition, as an 
application of Mercer’s theorem 


(5.3) 


Cor(/) = ^ Afe((^fe (g) (pk), 


where (Xk.fpk) are the eigenpairs of the associated Hilbert-Schmidt operator 


{D)BcP^ [ Coy{f){;y)cP{y)dyGL^{D). 
Jd 


Another way to obtain the decomposition (5.1) is a (possibly infinite) Cholesky decomposition of the two- 
point correlation function. Pivoting the Cholesky decomposition yields an extremely efficient approximation 

method, see e.g. [13 HU- 


5.2. Numerical calculation of a low-rank approximation of Cor(/). 


In general, the expansion (5.1) is infinite and has to be appropriately truncated for numerical computations. 
Let V := span{(pi : i = 1, 2,..., N} C L'^{D) be a suitable discretization of L^{D), e.g. a finite element space 
associated with a mesh of D. We are looking for a low-rank approximation of Cor(/) in the tensor product 
space V (S>V: 


(5.4) 


m X n \ / ^ \ 

Coiif){x,y) « ( E ^k,j>Pj{y)j e 1^0 K 

/c=l 2 = 1 ^ 7 = 1 ^ 


The unknown coefficient vectors in (5.4) can be computed as follows. Define the discrete correlation matrix 
C G as 

Oj = / / CoT{f){x,y)ipi{x)ipjiy)dxdy, i,j = l,...,N. 

J D J D 
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and the mass matrix G £ as 


Gi,j = / ipi{x)(pj{x)dx, i,j = 
JdD 


Then, it is easily seen that searching for a decomposition of the form (5.41 translates, in terms of matrices, 
into the search of an approximation Cm such that 


C Ki Cm = ^ 4^ with 4 = (4y)i=i,. 


.,N 


= G-4fc, 


fc=l 


is rigorously controlled (in a way yet to be defined). 


in such a way that the truncation error \\C — C, 

The best low-rank approximation in L^{D x D) is known to be the truncated spectral decomposition (5.3) 
(see e.g. |1H])- In the discrete setting, this corresponds to the spectral decomposition of G, which is a very 
demanding task. In particular, the decay of the eigenvalues {Afc} and thus the rank m to be reached for an 
accurate decomposition depend heavily on the smoothness of the underlying two-point correlation function 
Cor(/). Related decay rates have been proved in [28]. 

We suggest instead to employ the pivoted Cholesky decomposition in order to compute a discrete low- 
rank approximation of Cor(/), as originally proposed in [14] . It is a purely algebraic approach which is 
quite simple to implement. It produces a low-rank approximation of the matrix G for any given precision 
e > 0 where the approximation error is rigorously controlled in the trace norm. A rank-m approximation is 
computed in 0{m?n) operations. Exponential convergence rates in m hold under the assumption that the 
eigenvalues of G exhibit a sufficiently fast exponential decay, see El for details. Nevertheless, numerical 
experiments suggest that, in general, the pivoted Cholesky decomposition converges optimally in the sense 
that the rank m is uniformly bounded with respect to the truncation error e by the number of terms required 
for the spectral decomposition of Cor(/) to get the same error e. 


5.3. Low-rank approximation of the shape functional and its gradient. 


The basic idea is now to insert the state’s expansion (5.2) into the expressions of the expectation of the 


random shape functional and the associated shape gradient to derive computable expressions. In fact, it 
turns out that only standard solvers for boundary value problems need to be provided. Loosely speaking, 
this implies that, if one can compute the shape functional and its gradient for a deterministic right hand side, 
then one can also evaluate the expectation of the shape functional and its gradient for a random right-hand 
side. We illustrate this idea with the three examples introduced in Section]^ 

• Having decomposed Cor(tt) as Cor(M) = 'U'k®Uk, with each function Uk satisfying 


(5.5) 


-Auk = fk in D, 
Uk = 0 on dD, 


the mean value (4.2) of the Dirichlet energy can be computed as: 

D 


Moreover, its shape derivative is given by 


ye £ Bad, M\D)ie) = -^ [ 


duk 


On 


9 ■ n ds{x). 


• The mean value (4.51 of the L^-tracking type functional considered in Section 4.1.2 can be computed 
as 


M(D) = ^ f '^{uk- uq)^ dx 

B 7 . 
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with the Uk given by ( |5.5|) . As for the calculation of the shape gradient (4.6), we have to introduce 
the adjoint states pk G Hq(D), defined by 


J -Apfe = -{uk- uo) in D, 

I Pfe = 0 on dD, 

Thus, in view of the fact that Cot{pd,ud) = J2k Pk®Uk, we are led to the following formula for the 
shape gradient: 

V" e- (E If It) " ^ n *(.). 


• Last but not least, in the linear elasticity setting of Section |4^ expanding the correlation function 
of the surface loads Coi{g) = J2k9k ® 9k^ the two-point correlation Cor(u) satisfies the expansion 
Cor(u) = Uk ® Uk with Uk given by 

—div(Ae(Mfc)) =0 in D, 

Uk =0 on Tn, 

Ae{uk)n = gk on T^v, 

Ae{uk)n = 0 on T. 


Hence, the mean value Ai(D) of the compliance is given by 

M{D) = / Ae{uk) : e(uk) dx, 

k 

while its shape derivative reads 


V6» G 


0ad, M'{D){e) = Ae{uk) : 


9 ■ n dx. 


Remark 5. 


In the last example, one may be interested in the case that random body forces /(a;,w) are also 
applied to the system (which are, e.g. assumed to be uncorrelated with the surface loads g). Then, 
one has to perform two low-rank expansions Cor(/) = fk 0 fk and Cor(p) = 9i ® 9i which 
leads to an expansion of u of the form Cor(it) = J2k i '^k,i ® Uk,i with the Uk,i given by 

-div{Ae{uk,i)) = fk in D, 

Uk,i =0 on r^), 

Ae{uk,i)n = gi on Tat, 

Ae{ukj)n = 0 on T. 

The above formulae coincide with those for the multiple load objective functions (and their deriva¬ 
tives) proposed, e.g. in [3]. In contrast to this work, here, the different load cases are not known a 
priori, but originate from a low-rank approximation of the correlation function of the data. 
Hitherto, we have only been considering low-rank approximations of Cor(/) of the form Cor(p) « 
J2k9k ® 9k, where the gk are deterministic data functions, as they are naturally produced by the 
pivoted Cholesky decomposition. Notice however that the above discussion straightforwardly extends 
to the case of a low-rank decomposition of the kind Cor(p) « J2k i9k ^Hk Alfk® gk) (see the example 


of Section 6.2 for an application of this remark). 


6. Numerical examples 

We eventually propose two numerical examples which illustrate the main features of this article; both of 
them take place in the setting of linear elasticity as considered in Section |4.2| 
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6.1. Presentation of the numerical algorithm. 


When it comes to the numerical implementation of shape optimization algorithms, one main difficulty lies 
in the robust representation of shapes and their evolution. To achieve this, we rely on the level set method, 
initially introduced in |22j and brought to the context of shape optimization in mso]- 

The basic idea is to consider a shape D a.s the negative subdomain of an auxiliary ‘level set’ function 
(/) : > K, i.e. 

( (fiix) <0 if a; S H, 

\/x e < </>(a;) =0 if a: G dD, 

[ (j){x) >0 if a; G ^D. 

Thus, the motion of a domain D{t), t G [0,T], induced by a velocity field with normal amplitude V(t,x), 
translates in terms of an associated level set function •) as a Hamilton-Jacobi equation: 


( 6 . 1 ) 


^ + 'i/|V0|=o, tG (o,r), xgm^'. 


Hence, a (difficult) domain evolution problem is replaced by a (hopefully easier) PDE problem. Note that, 
in the present situation, V stems from the analytical formula for the shape derivative of the considered 
objective function Ai{D), which enjoys the structure (see Theorem]^: 


V0G0ad, M'iD){e)= Jvoe-nds, 


where 'Du is a scalar function. 

In numerical practice, the whole space is reduced to a large computational box Dq, equipped with a 
fixed triangular mesh T. Each shape D C Dq is represented by means of a level set function discretized 
at the vertices of T. In this context the elastic displacement ujo, solution to the linear elasticity system 


(4.10), which is involved in the computation of cannot be calculated exactly since no mesh of D is 


available. T herefo re, we employ the Ersatz material approach [1] to achieve this calculation approximately: 
the problem (4.10) is transferred to a problem on Dq by filling the void part Dq\D with a very soft material, 
whose Hooke’s law is eA with e ^ 1. 

All our finite element computations are performed within the FreeFem++ environment [24) . and we rely 
on algorithms from our previous works [SIIIQ], based on the method of characteristics, when it comes to 
redistancing level set functions or solving equation (6.1). 


6.2. Comparison between correlated and uncorrelated loads. 


This first example is aimed at appraising the influence of correlation between different sets of loads applied 
to the shapes. Let us consider the situation depicted on Figure (left): a bridge is clamped on its lower 
part, and two sets of loads Qa = (1, —1) and gt = (—1,1) are applied on its superior part, which may or may 
not be correlated. The actual loads are then modelled as a random field g{x,uj) of the form 

g{x,uj) = ^i{uj)ga{x) -i- ^ 2 {uj)gb{x), 
where the random variables and ^2 are normalized so that 

f ^,P(da;)=0, f e-P(dcc) = l, i = l,2. 

Jq Jq, 

The degree of correlation between ga and gt is measured by a := ^ 1^2 P(da;) where the case a = 0 

corresponds to uncorrelated loads. In this context, the correlation function Cor(g) G [T^(r 7 v)]‘^ x [L^(r 7 v)]‘^ 
naturally arises as a finite sum of pure tensor products (so that no low-rank approximation is necessary), 
and reads 

Cor(g) = ffa ® ffa + 5b ® + a (5a ® 5b + 56 ® ffa) ■ 

The mean value M.{D) of the compliance of shapes and its derivative can be calculated explicitly, along the 
lines of Section(see also Remark]^. 

We run several examples, associated to different values of the degree of correlation |a|< 1. In each situa¬ 
tion, an equality constraint Vol(Zl) = dx = 0.35 on the volume of shapes is enforced owing to a standard 
Augmented Lagrangian procedure (see [2T], §17.4). Starting from the initial shape of Figure(right), 250 
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iterations of the algorithm outlined in Section [6.1 [ are performed. The mesh T of the computational domain 
is composed of 12 141 vertices and 23 800 triangles; the CPU time for each example is approximately 12 
min on a MacBook Air with a 1.8 GHz Intel Core i5 with 4 GB of RAM. The resulting optimal shapes are 
represented in Figure]^ and the evolution of the objective function M.{D) and the volume Vol(Il) can be 
appraised on the histories of Figure [^ 

A huge difference in trends can be observed, depending on the degree of correlation between the loads 
(observe also the values of the objective function M{D) in Figure]^. Roughly speaking, as a gets closer 
and closer to —1, the shapes have to withstand the ‘worst-case’ of the situations entailed by ga and gb- 




Figure 2. Details of the test case of Section 6.2 (left) and initial shape (right). 



Figure 3. Optimal shapes obtained in the test case of Section |6)2l associated to degrees of 
correlation a = —1, —0.7,0,0.5, 0.8,1 (from left to right, top to bottom). 
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Figure 4. Convergence histories for the mean value (left) and the volume (right) in the 
test case of Section 16.21 


6.3. An example with a more complex correlation function. 


Let us turn to an example where the correlation function of the data is no longer trivial (i.e., it cannot be 
written as a finite sum of pure tensor products). The situation at stake is depicted on Figurea bridge is 
clamped on its lower part and (random) surface loads g = ( 31 ,( 72 ) are applied on its top. 

We study three different scenarii, corresponding to surface loads 3 * = ( 31 , 32 ), i = 1,2,3. For the sake 
of simplicity, in all three cases, the horizontal and vertical components g\ and g^ are uncorrelated; the 
associated correlation functions are given by 


Va;,3 G Tn , 


Cor{gl){x,y) = 10 ^ h+ ^ ^ 

Cor(3^)(a;,3) = 10® fc+ ^ 


where the superscript stands for the positive part, the characteristic length I is taken as I = 0.1, and the 
functions hi and ki {i = 1,2,3) are defined as (see also the graphs in Figure]^ 


/ii(t) = 1-4 0 , 

fci(t) 

h2{t) = 2t{l — t) + -, 

^ 2 ( 1 ) 

hsit) = 1, 

ksit) 


16(1-1)^ 

, if t < 

1 

2 ’ 

16(1-!)^ 

, else. 



if t < 

24(t-|)(t-i), 

else. 

24(t-|)(i-i), 

if t < 

24(1-1) (1-1), 

else. 


1 

2 ’ 


1 

2 ’ 


Loosely speaking, these correlation functions show a decreasing dependence on the distance \x — y\ between 
two points x,y G F^v, and the factors hi,ki, which depend only on the average position {x + 3)/2, mimic a 
variable intensity of the loads according to the spatial location. Note that the pivoted Cholesky decomposi¬ 
tion, as described in Section is used to obtain low-rank approximations of these correlation functions of 
the form (5.1), and we retain the first five terms in each expansion. 

The objective function of interest is, again, the mean value A4(D) of the compliance of the structure. A 
constraint Vol(Zl) = 0.75 is imposed on the volume of shapes and 250 iterations of the algorithm outlined 
in Section |6.1| are performed in each situation on a computational mesh composed of 5 752 vertices and 
11202 triangles which requires a CPU time of approximately 15 min. The resulting optimal shapes and 
convergence histories are reported on Figures and respectively, showing very different trends depending 
on the particular situation. 
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Figure 5. Graphs of the functions hi (left) and graphs of the functions ki (right). 



Figure 6. The setup of the test case of Section 6.3 (left) and initial shape (right). 
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